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Abstract. 

A new lattice method is presented in order to efficiently solve the electrokinetic 
equations, which describe the structure and dynamics of the charge cloud and the flow 
field surrounding a single charged colloidal sphere, or a fixed array of such objects. 
We focus on calculating the electrophoretic mobility in the limit of small driving field, 
and systematically linearise the equations with respect to the latter. This gives rise to 
several subproblems, each of which is solved by a specialised numerical algorithm. For 
the total problem we combine these solvers in an iterative procedure. Applying this 
method, we study the effect of the screening mechanism (salt screening vs. counterion 
screening) on the electrophoretic mobility, and find a weak non-trivial dependence, as 
expected from scaling theory. Furthermore, we find that the orientation of the charge 
cloud (i. e. its dipole moment) depends on the value of the colloid charge, as a result 
of a competition between electrostatic and hydrodynamic effects. 



PACS numbers: 47.ll.-j, 47.57.jd, 47.57.E-, 47.57.J-, 47.57.-s, 82.45.-h, 82.70.Dd, 
82.70.-y, 05.60.Cd 
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1. Introduction 

The interplay between electrostatic and hydro dynamic interactions is of high importance 
for the understanding of a wide range of biological, chemical and physical systems, since 
in almost all situations where a solid is brought in contact with a liquid, a difference 
in the electric potential occurs due to association or dissociation of charges or the 
orientation of molecules at the surface. 

Charged solid colloidal spheres in suspension in an aqueous solution containing 
counterions and salt ions will be surrounded by a cloud of oppositely charged ions. 
This cloud, typically called the electric double layer, is responsible for screening the 
electrostatic potential. If an external electric field is acting on the system, the charged 
spheres start to migrate in the direction of the oppositely charged electrode and 
the surrounding cloud will be deformed and becomes anisotropic due to the electric 
field and also to the friction between the ions and the fluid. This phenomenon is 
called electrophoresis and the corresponding transport coefficient is the electrophoretic 
mobility determined by the balance of electric driving and hydrodynamic frictional 
force acting on the sphere. It is defined as the proportionality constant between the 
constant velocity u of the particle and the external driving field E^xt in the linear 
regime, i. e. for small driving fields, 

U = flEext ■ (1) 

Efforts have been made to study electrophoresis by experimental methods [H El El HI E] 
as well as by analytical and numerical calculations over the decades [HI [71 El [9l [TOl [TTl 
[la [HI [11 [151 [161 [HI [TBI [191 EQl Ell Due to the complicated many-body 

nature of the problem a comprehensive quantitative theoretical understanding is still 
lacking. 

An important milestone was the theoretical investigation by O'Brien and White 
|13] . later known as "standard electrokinetic model". Starting from the dynamic Mean- 
Field equations that describe the interplay between the convection-diffusion dynamics 
of the ion clouds, the solvent fiow field, and the electrostatic forces in the system, they 
studied a single charged colloidal sphere in infinite space, subject to an infinitesimally 
weak homogeneous external electric field, with respect to which the equations are 
linearised. This problem exhibits full spherical symmetry, and hence its numerical 
treatment reduces to the solution of a one-dimensional ordinary differential equation, 
which finally allows the calculation of fi. 

This approach has a fairly broad but not unlimited range of applicability, whose 
conditions may be summarised as follows: Firstly, the Mean-Field theory as such must 
be justified, and this implies weak ion-ion correlations, which is typically the case for 
single-valence ions at room temperature. Furthermore, treating the problem within the 
framework of a single-colloid theory requires that the ionic clouds essentially do not 
overlap, and all non-trivial values of charge density, electrostatic potential, and fiow 
velocity are confined to the double layer, as a result of electrostatic and hydrodynamic 
screening. This means that the salt concentration has to be fairly large, and the 
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electrostatic screening is dominated by the salt contribution: Note that for a single 
sphere in infinite space the counterions have all entropically "evaporated", and hence 
must be ignored in the theory. 

It is exactly this latter condition that is violated in recent experiments on colloidal 
electrophoresis [H |22], which have deliberately focused on the limit of low salt. In 
Ref. |22] it was shown that some (certainly not all) experimental observations in 
that regime can be explained by assuming that the effects of finite colloid volume 
fraction (and corresponding finite counterion concentration) can be modeled by simply 
studying a single colloidal particle in a finite simulation box with periodic boundary 
conditions, which automatically gives rise to a finite colloid volume fraction and the 
correct corresponding counterion concentration. This investigation was done by studying 
a system of charges by Molecular Dynamics, while the solvent hydrodynamics was taken 
into account by a Lattice Boltzmann background. However, the computational effort of 
such studies turned out to be large, such that it was neither possible to obtain highly 
accurate results, nor to vary parameters systematically over a broad range. 

For these reasons, we develop a new approach in the present paper: On the one 
hand, we wish to study precisely the same physical situation as in Ref. [22], i. e. a 
single colloidal particle in a finite box with a finite counterion concentration, possibly 
with added salt, while on the other hand taking full advantage of the O'Brien and 
White approach, which means that we study the Mean-Field equations, combined with 
linearisation with respect to the driving field. In summary, our work is nothing more 
and nothing less than the finite-volume generalisation of Ref. [I3], whose results are 
directly recovered in the special limit of large salt concentration, where the double layer 
is significantly smaller than the box. Our method is based upon a full three-dimensional 
calculation, where the partial differential equations are discretised on a lattice. Since 
these equations are mathematically fairly different, we develop specialised solvers for 
each equation. The colloidal particle acts essentially as a boundary condition, and 
hence it is clear that the method can be easily generalised to a fixed array of particles 
or any other object that is periodically repeated in space — strictly spoken, we are 
studying a colloidal crystal, as a result of the periodic boundary conditions of the box. 
Within this setup one may vary several important parameters like volume fraction and 
averaged ion concentrations, and study the behaviour of fi fairly accurately. 

Returning to the issue of the limitations of the present approach, we notice that 
the present model differs from the Molecular Dynamics / Lattice Boltzmann (MD/LB) 
model not only with respect to its computational cost, but also in terms of the modeling 
of the finite size of the ions. The present model clearly assumes point ions and neglects 
any ion-ion correlations, which, however, due to packing effects, are important if the 
colloid size is not very large compared to the ion size [151 [IH] • On the other hand, these 
correlations are fully present in the MD /LB model. However, it should be noted that the 
MD/LB model is, for computational reasons, very limited in terms of the ratio of colloid 
size to ion size, which cannot have values much in excess of 10. Therefore, the MD/LB 
model will probably, in comparison to experiment, overestimate the effects of packing 
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and ion correlations, except for quite small colloids. We therefore study here the other 
extreme, which should be reasonable for large colloids. While there are attempts to 
include the finite size of ions into more generalised Mean Field theories, both for statics 
[I5| [T6| 125] and dynamics [151 EH], it is not immediately obvious how to incorporate 
these formalisms into a three-dimensional code that is strictly confined to a finite box 
with a well-defined and conserved number of ions. Furthermore, the present model 
should be viewed as just a first step to the development of a more general simulation 
program that is able to study multi-colloid systems. At finite volume fraction, colloid- 
colloid correlations are quite important, and they are of course not taken into account 
by our "cell" or "colloidal crystal" model. For the statics, the importance of such 
correlations has been pointed out, e. g., in Refs. [26| [27] : in the dynamic case one 
expects a probably even stronger effect, since here not only the electrostatic interactions 
are insufficiently screened, but the hydrodynamic interactions as well. A multi-colloid 
simulation model with explicit ions is however very likely to be computationally too 
expensive, except for very moderate scale separations between macro- and micro-ions, 
both with respect to length and with respect to charge. The present work is intended 
as a first step in an attempt to overcome at least the former of these limitations, by 
assuming full scale separation between length scales at the outset. Of course, numerically 
this scale separation is anything but perfect, due to limited grid resolution; however, 
discretised field theories seem to be more amenable to extrapolation procedures than 
particle models. 

In Sec. [2] the theoretical model is introduced, including the Mean-Field approach 
and the linearisation of the equations. Furthermore, the equations are reformulated in 
dimensionless units. The computational method is briefly discussed in Sec. [3], where we 
describe the iterative combination of the specialised solvers and discuss the particular 
choices for the numerical methods. In Sec. H] some interesting results from this method 
are presented: The dependence of the electrophoretic mobility on the details of the 
screening mechanism is analysed, and the reversal of the field-induced dipole moment 
of the ion cloud surrounding a weakly charged colloid is elucidated. Section [5] concludes 
with a brief summary. 

2. Theory 

2.1. Electrokinetic equations 

Electrophoresis is a result of the balance between electrostatics and hydrodynamics. 
Within a Mean-Field picture the system is described in terms of ion concentration fields 
Cj, electrostatic potential ip and the fiow velocity field v. Cross-correlations between 
salt ions as well as thermal fiuctuations are neglected. 

The Poisson equation couples the concentration fields to the electrostatic potential. 




(2) 
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Here, e is the dielectric constant, e denotes the elementary charge and Zi is the valence 
of the ionic species, where the subscript i indicates the different ionic species in the 
system. Counterions, which assure the charge neutrality of the system, are denoted by 
the index i = 0. The charged colloids are taken into account via boundary conditions. 

The dynamics of the concentration field q is described by a continuity equation, 
where the total current density is a combination of a diffusive term, a convective current 
and the current resulting from the electric force. One thus obtains a convection-diffusion 
equation, known as Nernst-Planck equation, 

dtCi = V ■ ^AVci + -^^ezi(yij)ci - vci^ . (3) 

Here, Di is the diffusion constant of the ionic species i and ksT denotes the thermal 
energy. The ion mobility is given by Di/iksT) due to the Einstein relation. 

Electric forces and viscous forces are balanced in the Stokes equation which 
describes zero Reynolds number incompressible hydrodynamics, 

V ■ t; = , (4) 

pdtv = -Vp + r]V'^v - e{Vij)^ZiCi, (5) 

i 

where p is the mass density of the fluid, p its pressure field and rj the fluid viscosity. 
In the stationary state, the system of equations is thus summarised as [28] 

= VV' + -e XI , (6) 

i 

= V ■ (^AVq + -^eziiy^Yi - vc}j , (7) 
= - Vp + rjV'^v - e{V^) ^ ZiCi , (8) 

i 

Q = V -v. (9) 

Notice that the stationary formulation is not manifestly Galilei invariant, but rather 
selects one particular frame of reference (the rest frame of the colloidal particle) in which 
it is valid. If a colloidal sphere would move relative to the chosen inertial frame, the 
local ionic concentration would change with time, and hence a stationary solution would 
not exist. It is this restriction which either confines the method to single-colloid studies, 
or forces us to impose a somewhat unphysical "rigid-body" constraint between the set 
of colloidal particles. The mobility however is measured in the system's center-of-mass 
reference frame. In other words, the center-of-mass velocity of the system must be 
taken as the velocity that determines p. As a matter of fact, it turned out that it is 
most convenient to solve the Nernst-Planck equation in the colloid rest frame, while the 
Stokes equation is best solved in the rest frame of the center of mass. Therefore, one 
always needs a trivial Galilei transform when switching from one equation to the other. 
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2.2. Dimensionless formulation 

An important length scale in the theory of charged systems is the Bjerrum length, which 
results from the balance between electrostatic and thermal energy: 

The Stokes mobility of a sphere of radius Ib and elementary charge e provides a 
natural unit for the electrophoretic mobility: 

bTirjlB 

and the dimensionless reduced electrophoretic mobility is defined as 

/^red = — • (12) 

The natural energy scale is the thermal energy /c^T, and together with the elementary 
charge e this yields a dimensionless electrostatic potential 

^ = ^^. (13) 

Introducing a second length scale (see below), such that the gradient is rescaled via 

V = -V, (14) 

K 

the electrokinetic equations are nondimensionalised as 

Q = V^i, + Y.ZiCi, (15) 

i 

= V ■ [DiWc, + b,Zi{y^)ci - vc^ , (16) 
= - Vp + '^-V^v - (V^) ^ z,c, , (17) 

i 

= V-i;, (18) 

where the dimensionless parameters and variables are summarised in Tab. [TJ 

Note that the transformation from physical to reduced units, as outlined in Tab. [H 
is valid for any arbitrary choice of the length scale . However, a physically motivated 
choice results from the finite-volume version of linearised Poisson-Boltzmann theory 
(Debye-Hiickel theory). We therefore choose k to be the Debye screening parameter, 

t^ = AnlBY,^'T7^ (19) 



where all ionic species (including the counterions) contribute, V is the volume of the 
system (actually the volume that is available to the ions, i. e. box volume minus colloid 
volume), and Ni the number of ions of species i. For simplicity the tilde will from now 
on be omitted, with the understanding that all parameters are given in reduced units. 
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parameter 


physical units 


dimensionless formulation 


Bjerrum length 


2 

I — 




screening parameter 


K^ = A7rlsE^^ff 




electrophoretic mobility 


fJ' 


t^red e ^ 


spatial position 


r 


r = KV 


spatial derivative 


V 


v = ^v 


electrostatic potential 


i; 


^ = k^^ 


electric field 


E 


TP ^ TP 


ion concentration 


Ci 


7. — ^-^Ib^. 


colloid charge 


Ze 


Z = A-kIbuZ 


number of ions 




Ni = AnlsKNi 


fiow velocity 


V 


V = ^I'«> 

KKBl 


pressure 


P 


~ _ AttIb „ 
P - K^kgTP 


diffusion constant 


Di 


f\ _ GttIbv n 

— UbT 



Table 1. Summary of all parameters in physical units and their reduced counterparts. 
2.3. Linearisation 

The high nonlinearity of the Mean-Field equations causes two problems. Firstly, 
it is difficult and computationally expensive to solve a coupled system of nonlinear 
differential equations. Furthermore, the electrophoretic mobility is well-defined (i. e. 
independent of the driving field) only in the linear regime. Consequently, if a fully 
nonlinear solution of the equations is obtained, an extrapolation to zero driving field 
is required. The second problem can be avoided completely, and the first one at least 
reduced, by a linearisation of the equations in terms of the driving field [I3l |28]. This 
can be done by a formal expansion with respect to a small parameter e, corresponding 
to the strength of the external field. All fields in the system have a regular expansion 
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in e, and hence may be written as 

c, = cf ) + ecf ) + 0{e') , (20) 
^ = ^(0) + eV^(i) + (9(^2) ^ (21) 

V = ev(i) + C(e2) , (22) 
p = p(o) + ep(i) + 0(e2) . (23) 

note that for e = 0, i. e. in the absence of external driving, the system is at rest, such that 
the zeroth-order velocity vanishes. We now insert the expansion into the electrokinetic 
equations; noticing that for e = all ionic currents must vanish, one obtains 

• in zeroth order perturbation theory: 

= vV^ + ^z,cf\ (24) 

i 

= v(z.^(°)+lncf^) , (25) 
= - Vp(°) - (VV^(°)) ^i^T ■ (26) 

i 

• and in first order: 

= V-{AVcf)+A;2.(V^A«)cf 

+ A;^.(V^(°))c«-««cf I, (27) 



= - Vp(i) + \w^v^^^ - (VV^(i)) ^^,cf - (VV^(°)) J]z,cf) , (28) 

I I 

= V ■ v^^^ , (29) 

Q = sj^^m + Y^z4\ (30) 

i 

The zeroth order only contains the electrostatic potential of the unperturbed ion 
clouds; hence this order is identical to standard nonlinear Poisson-Boltzmann theory. 
Equation [26] is just an equation to determine the zeroth-order pressure, which is of no 
interest to us; it can therefore be simply ignored. The first order consists of a coupled 
set of linear equations; hence the only nonlinearity that remains is the equilibrium 
Poisson-Boltzmann problem, which is simpler than studying the original full set of 
nonlinear dynamic equations. 

In the first-order equations, the external field is taken into account by decomposing 
the potential ■z/^^^-* into a periodic part and one part corresponding to the constant electric 
field 

^(1) = ^'(1) + , (31) 



such that 



VV'« = -Y.'^<^^ (32) 

i 

V'V"^'^ = , (33) 

V^"(i) = -E,,t. (34) 
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Hence one may write the first-order equations more explicitly as 
= V ■ { AVcf ^ + DiZ,{V^'^'^)cf^ - DiZ,E,,tcT 

+ A^^(V^A(°))cf)-vWcf }, (35) 

= _ VpW + 2^2^(1) _ ^^^,(1)^ v-^^^(o) 
o 

i 

+ 5^ z,cf)-(V^W)5;]z.cl^ (36) 

i i 

= V ■ 'y(^) , (37) 
= vV(^) + ^z,cS^^ (38) 

i 

It should be noted that, as a result of the perturbation expansion, the first-order 
convection-diffusion equation now contains sources and sinks (the terms proportional 
to cf''). However, since these terms all have the form of a divergence, there is neither a 
total flux of matter into the system, nor out of it, as it should be, since mass conservation 
must hold at each order of the expansion separately. 

The reduced electrophoretic mobility is then finally calculated as 

l^red = 1 , (39) 

\-^ext\ 

where w*^^^ is the constant velocity of the colloid in the system's center-of-mass reference 
frame, i. e. (assuming no-slip boundary conditions) the value of the flow velocity field 
at the surface of the colloidal sphere, 

1,(1) =-yW(i?), (40) 

with R the radius of the particle. This mobility is strictly independent of the strength 
of the external driving field. 



3. Computational method 

3.1. Iterative procedure 

The linearisation of the problem divides the challenge of solving the electrokinetic 
equations into two different subproblems. For the zeroth order, a solution of the fully 
nonlinear Poisson-Boltzmann equation must be found. The first order consists of a set 
of linear equations, where the zeroth order fields only occur as prefactors. Here we need 
to solve the Poisson equation, the convection-diffusion equation and the incompressible 
Stokes equation. Each particular problem can then be dealt with by using a specialised 
solver, which is some finite-difference scheme on a regular lattice with periodic boundary 
conditions; details for each solver will be outlined below. Finally all methods are 
combined via iterative loops, as sketched in Fig. [H 
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Input 1: 

Equilibrium 
parameters 




YES 



Nernst-Planck 

(one time step) 




NO 



Figure 1. Schematic illustration of the iterative algorithm for solving the 
electrokinetic equations. 



It turned out that the convergence of the method is improved by not using the full 
new velocity field for the next iteration, but rather a convex combination of the result 
of the previous iteration and the most recent result of the Stokes solver v*: 



(1) 



< w < 1 ; (41) 

in practice we used u = 1/2. The iteration procedure yields a sequence of reduced 
mobilities /x^^^, z = 1, 2, . . .; the iteration was terminated as soon as the relative residual 



(i) 
f'red 



(i-l) 
f'red 



(i) 
l^red 



(42) 
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dropped below the value 10 ^. 
3.2. Poisson-Boltzmann equation 

For the zeroth order a solution of the fully nonlinear Poisson-Boltzmann equation is 
required. Here, a simple and unconditionally stable lattice algorithm has been used, 
based on a constrained variational approach of the Poisson-Boltzmann equation. This 
method has been discussed in detail in Ref. [29] and is hence only be briefly summarised 
here. It should be noted that this method has prompted other groups to develop similar 
ideas even further [25]; however, such more recent implementations have not been used 
here. 

Following the ideas of Maggs and Rosetto [30] and re-formulating the equations 
in terms of the electric field instead of the electrostatic potential, Eqs. [2l] and [25] are 
written as 

V-E = z,cu (43) 

i 

V X E = 0, (44) 
ZiE =VlnCi. (45) 

These equations are recovered as the Euler-Lagrange equations of a constrained free 
energy functional of the form 



T= / fdV, (46) 
Jv 

f =^E' + J2c^lnc,-^i\/■E-J2^^c,)-J2^^^^ic^-y), (47) 

i i i 

where the electrostatic potential ip and the chemical potential /ij of ionic species i 
occur as Lagrange multipliers taking into account Gauss' law and mass conservation, 
respectively. Ni is the total number of ions of the species i and V denotes the system's 
volume. Note that the formulation in terms of the electric field assures that the solution 
of the Poisson-Boltzmann equation is a true minimum of the functional. Applying a 
Yee discretisation [31], i. e. associating scalar fields with the sites, polar vectors with the 
links, and axial vectors with the plaquettes of a simple-cubic lattice, the functional can 
be minimised making local moves between adjacent nodes along a link and local field 
updates on the plaquettes. If the system has been initialised such that the constraints 
are fulfilled, those local moves never leave the constraint surface. Moreover, the update 
rules can be optimised, such that the functional value is decreased in every iterative 
step, and the method will run ultimately into the one and only minimum. For further 
details the reader is referred to Ref. [29] . 

3.3. Poisson equation 



The Poisson equation for a given charge density can be solved efficiently by Fast Fourier 
Transform. We expand the potential and the charge density in terms of Fourier series 
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via 



with 



ijj{r) = ^ '4){k) exp[-ik ■ r] , (48) 

k 

P(^) = P^^") ^M-ik ■ r] , (49) 

k 

k = 27r{k/L^,l/Ly,m/Lz) , kJ,meZ. (50) 



^ = -nP- (51) 



Here, Lr X Ly X is the dimension of the computational domain. The solution of the 
Poisson equation in Fourier space is then given by 
1 
k 

For consistency reasons we use a discretised version, i. e. a lattice Green's function, 
instead of the continuum Green's function [32]; the discretised counterpart to Eq. [5T] 
reads 

ip^k, I, m) = -^—p{k, I, m) (52) 



with 



eL,m = ^ <! 3 - cos ( 27r-^ ) - cos ( 27r-^^ - cos (2tx^ ) \ , (53) 



NJ V Nyj V iV. 



where a denotes the lattice spacing and A^^ = Lx/a etc. Back-transformation finally 
yields the desired electrostatic potential in real space. 



3.4- Stokes equation 

The stationary incompressible Stokes equation has the form 

V-v{r) =0, (54) 

-Vp(r) + r,V^v{r) = - f^,,{r) , (55) 

where f^^t is an external force density and r] denotes the fluid viscosity (which takes 
the value 2/3 in our reduced unit system). For the purposes of the present paper, one 
should view f^^^ as the force density generated by the electric field and the charges. 
More precisely, we include in f^^t electric forces that come from the ion clouds, but 
also the force density that is generated from the fixed charges of the immersed body (or 
bodies). Since the total system is charge neutral, the total force on the system vanishes, 
even in the presence of external driving: 

/ rfV/,,,(r)=0. (56) 
Jv 

Again, it is obvious that this statement holds for the full nonlinear theory, and this 
implies that it must also hold separately at each order of the perturbation expansion. 
We now assume that the immersed bodies are not moving relatively to each other (cf. 
the remark at the end of Sec. 12. ip . and also not rotating. Under these circumstances. 
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the flow field can be calculated rather straightforwardly, making use of the idea of 
replacing the differential equations by equivalent integral equations [33], [M]. This is 
done as follows: The boundary condition for the fluid is given by a unique constant 
(but unknown) velocity on the surface. To assure this boundary condition, one may 
introduce an artificial "reaction force density" f^eac located on the surface (with units: 
force per area). This force density needs to be determined self-consistently such that the 
superposition of the flow fields generated by the external forces and these reaction forces 
satisfies the boundary condition. This is in spirit quite analogous to the electrostatic 
problem of a metallic surface, where the problem of finding a constant electrostatic 
potential at the surface is solved by determining an appropriate induced charge density. 
It is clear that this reaction force density cannot exert a net force onto the system, and 
hence we know 

'"rffi/,_(r) = 0; (57) 

here Q denotes the surface. 

We thus can write the total flow field i; as a superposition of vi, the contribution 
from the external force density, and V2 coming from the reaction force density: 

v{r) = vi{r) + V2{r) (58) 

v,{r)= f d\' T{r-r')f,,,{r'), (59) 
Jv 

v,{r) = [dQ'T{r- r') /_(rO , (60) 



where T is the Green's function of the Stokes equation. 

For an infinite fluid this Green's function is well known and given by the Oseen 
tensor (see e. g. |35]). In Fourier space it is given by 

where / denotes the unit tensor. In a finite box, the same form still applies; however, 
one needs to take into account that only wave vectors k occur that are compatible with 
the box periodicity. Furthermore, fc = must be excluded, since the problem is solved in 
the system's center-of-mass reference frame. This yields for the real-space counterpart 

here V denotes the total volume of the box. It should be noted that we can thus consider 
Vi as known, while we do not know V2 and /,,eac- 

Now, picking two points and r^e/ that are both located on the surface, we know 
that their velocity must be identical: 

Vl{r^) + V2{rn) = ViiVref) + V2{rref), (63) 

V2{rn) - V2{rref) = vi{rref) - vi{rn), (64) 
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I dn' (t {m - rn')- T {Tref - rn')) freadrn') 

= vi{rref) - vi{rn). (65) 

In a discretised version the Q' integral is replaced by a sum. If we now view r^e/ as an 
arbitrary but fixed reference point, while we vary vq, we may view Eq. |65]as a system 
of linear equations to determine freac- If the number of surface points is M, then the 
number of equations is 3M, while the number of unknowns is 3M as well. However, 
three of these equations are redundant, since at the point = Vref only the trivial 
information = is obtained. Instead, we need to use Eq. [57] as last set of equations to 
obtain a unique solution. In practice, the set of equations was solved numerically using 
the standard BiCGStab procedure [36| [37]. 

For discretisation, we again use a simple-cubic lattice with spacing a. For 
consistency reasons, we need to use the discrete version of the Oseen tensor, analogously 
to the lattice Green's function of the Poisson equation (see Eqs. [52] and [53] and Ref. |32j). 
The discretised Oseen tensor is derived by applying a midstep finite-difference scheme in 
real space, and doing the corresponding discrete Fourier transform with integer indexes 
k, I, m: 

T{k, /, m) = ( J (gg) 

V<.k,l,m \ ^k,l,m J 

with 

"M,™ = ^ f ta (2.^) , Bin (2xi-) . sin (2.^) | . (67) 
= I {3 - cos (2.±) - cos (2.±) - cos (2.^) } . (68) 

3.5. Convection- diffusion equation 

Equation [3S] is the stationary limit of a convection-diffusion equation of the general form 

— + V ■ w{r) \ c{r, t) = DV'^c{r, t) + S{r, t) . (69) 

Here c(r, t) denotes the ionic concentration field of first order at the spatial position 
r and time t. The convective term is not the fluid velocity, but rather the coupling 
term of the first order ion concentration with the zeroth order electrostatic potential, 
w{r) = — D ip^'^^r) . The source term S{r,t) contains all other terms independent of 
the first-order ion concentration field. Equation [35] is a conservation law, and therefore 



d^rc{r,t) = const. , (70) 

J d^rSir,t) =0. (71) 

The latter equation states that no ionic particles are produced or annihilated during 
the process. The discrete counterpart of such an equation is a Master equation, the 
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coefficients of which need to be adjusted such that its continuum hmit recovers Eq. EHl 
A detailed and systematic derivation of such an algorithm shall be published separately 
[38] . In the present work, we use the simplest version, which is a nearest-neighbour 
model on the simple-cubic lattice, which has lowest-order accuracy, and present the 
algorithm without proof. The concentration fields are initialised as zero at every lattice 
site and then propagated using a Master equation of the form 

c(r, t) = Y^ Ai{r - aAi) c{r - aA^, t - h) + hS{r, t ~ h) , (72) 

i 

where h is the discretisation time step, r denotes the lattice sites and a is the lattice 
spacing. Aj denotes a dimensionless lattice vector connecting r with one of its 
neighbours, such that r + aAj is again a lattice site. For our simple three-dimensional 
nearest-neighbour model, i = 1, . . . ,6, and the transfer coefficients are given by 

A{r) = ^ (l + ^^(r) ■ A.) , (73) 
while the diffusion constant is 

This means that the continuum limit is obtained by Taylor expansion up to second order 
with respect to only one variable a, since the time step is not a variable that can be 
picked independently from a, but rather is given hj h = a? /{QD). 

4. Numerical results 

4.1. Parameters 

From now on, we will return to the notation of Sec. 12.21 and, for reasons of clarity, 
distinguish between dimensional and reduced quantities. 

The electrophoretic mobility ^red is a dimensionless quantity, and hence it can only 
depend on dimensionless parameters as well. In Ref. [23] we discussed, within the 
framework of the present Mean-Field treatment, a finite system with added salt, where 
all ion types have the same properties, i. e. all ions are monovalent and have all the 
same friction coefficient. We then found as one possible set of dimensionless parameters: 

(i) the reduced charge 

Z = Z^ = ^, (75) 

(ii) the rescaled colloid radius R = kR, (iii) the rescaled diffusion constant of the ions 
D (as defined in Tab. [1]) and (iv) a dimensionless quantity /o that specifies the fraction 
of counterions (species zero) relative to the salt ions. In general fi is the fraction of the 
ionic species i relative to all ions in the system. 
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and it is easily shown that /, is nothing but the volume-averaged concentration of ionic 
species i, in the reduced unit system of Tab. [1] For a system with only one salt species 
that has only monovalent ions, we know /i = /2 (these fractions refer to the salt ions), 
and the sum rule Y^- zffi = 1 implies that only one non-trivial parameter /g is left. 

As already mentioned in Sec. 12.21 a definition of k that is fully consistent with the 
finite-volume version of linearized Poisson-Boltzmann theory requires that the volume 
V is defined as the volume available to the ions. For a box of dimension L x L x L and 
colloidal spheres of radius R this means 

471 

In our notation, is the number of colloids and Z is their charge, while A'o and zq are 
the number and valence of counterions, respectively, such that charge neutrality implies 

NZ = -NqZo. (78) 

The colloid volume fraction is thus given by 

or 

A-kNR^ $ , , 

T— = T^- 

Therefore the correct relation between volume fraction $ and /o differs from the 
expression in Ref. [23] by a small correction term. Inserting Eqs. [191 ESI ES] and 
1751 one finds after a few lines of algebra 

^ ^0/0/..nN2 _ ^0/052 /oiN 



1 - * 3Z 3Z 

Furthermore, a dimensionless resolution d is defined such that for given d a sphere 
is always discretised by the same number of lattice sites: 

d=|, (82) 

where d denotes the lattice spacing in reduced units and R is the radius of the particle, 
also in reduced units. 



4-2. Comparison with previous results 

Figure [2] studies the reduced mobility for a system consisting of one colloidal sphere in a 
box, where the volume fraction $ ~ 7.07- 10~^ as well as the reduced charge Z = 8, 10 are 
kept fixed. The counterions are monovalent. R is varied by changing /o (cf. Eq. [HT]) . i. e., 
by varying the salt content (one species, all salt ions monovalent). The calculations were 
done with a reduced diffusion constant = 1.5 for all ionic species, and the resolution 
was kept fixed at the value d = 0.07. One clearly sees that the mobility systematically 
decreases with R, which is easily explained by the corresponding increase of electrostatic 
screening. Note that the present representation is given for constant Z, which differs 
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Figure 2. Reduced mobility as a function of R. The open symbols are results from 
Ref. [22|. 



from the classical calculations [HI HDl [HI [I2l [13] that keep the zeta potential fixed, i. e. 
the electrostatic potential at the colloidal surface. 

Furthermore, Fig. [2] shows also data from Ref. [22] (open symbols). The circles 
and squares are simulation results obtained from the Molecular Dynamics / lattice 
Boltzmann (MD/LB) raspberry model [39]. In the simulations a single colloid of charge 
Z = 20 {Z = 30) and radius Rq = 3 {Rc = 5) is surrounded by counterions. Both 
systems were studied with a Bjerrum length of Ib = 1.3, resulting in a reduced charge 
oi Z = 8.5, comparable to our value. The triangles are experimental results for latex 
crystals in a bcc structure with a particle size of R = 34nm in a deionised aqueous 
suspension. The effective charge is quoted as Z ~ 450 determined from conductivity 
measurements [IQ], resulting in a reduced charge of order Z ~ 10. In all cases of Ref. 
[22] the colloidal size and charge were fixed and salt ions were absent (except for the self- 
dissociation of water). The screening parameter, or R, was hence changed by varying $, 
while /o was kept constant (cf. again Eq. [HT]) . This is qualitatively different from our 
numerical calculations, where rather $ is kept constant and /o is varied. Nevertheless, 
the results seem to agree quite nicely, and within the accuracy of the data it seems that 
it does not matter whether the screening is salt-dominated or counterion-dominated. 
In the following subsection, we will put this question under more detailed and more 
accurate scrutiny. 

Besides the question of the screening mechanism, there are also additional 
differences between our calculation and Ref. [22], in view of which the observed small 
discrepancies are hardly surprising: The MD/LB simulations use a slightly different 
reduced charge, and also the reduced diffusion constant of the ions is probably somewhat 




Figure 3. Reduced mobility as a function of the dimcnsionless (reduced) zeta potential 
(determined from the solution of the Poisson-Boltzmann equation) for R = 8, d ~ 0.07, 
D = 1 and fo — 0.01. The solid line is the result taken from Ref. [T5] . 



different. The inffuence of the diffusion constants on /i has so far not been thoroughly 
studied; in the next subsection it will be shown that /i increases with the Di, as one 
might expect. For the experiments, the situation is even less clear, since there is a 
considerable amount of ambiguity in the mapping of the effective charge of the real 
physical system to the bare charge in the Mean-Field calculations. 

Finally, Fig. [2] presents also a comparison with two theoretical results. Firstly, the 
Smoluchowski limit [6] of the reduced mobility is given by 

s 3 

/^reT = , (83) 

where (red is the reduced (dimensionless, cf. Tab. [1]) zeta potential, i. e. the electrostatic 
potential at the colloid surface. This can be easily calculated within our approach; it 
is just a result of our Poisson-Boltzmann solver, where we take for C the potential 
difference between colloid surface and box boundary. Secondly, within an approximate 
numerical theory for a single colloid in an infinite salt solution, Wiersema et al. [11] 
have tabulated values for fired as a function of (red and R, while the influence of Di is 
stated there to be fairly small. We can therefore use our values for R and (red to also 
compare with that theory, using linear interpolation. While the Smoluchowski limit is 
reached for values R ^ 40, the data obtained from the work of Wiersema et al. describe 
our numerical results quite well over the full range. 

Since our method is the extension of O'Brien and White's work [13] to systems with 
finite volume fraction and a finite amount of counterions, it should produce identical 
results in the limit of strongly salt-dominated screening (/o — ?■ 0), where the ionic 
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cloud of the colloid does not overlap with those of its periodic images. A quantitative 
comparison is however hampered by the fact that Ref. [13] does not tabulate its results, 
but only provides plots, and, more importantly, that the values of the ionic diffusion 
constants are not quoted. It seems however that D = 1 for all ionic species is reasonable; 
at least we do find quite good agreement between our calculations and Ref. [I3] for this 
value, as demonstrated in Fig. |3l 

4-3. Salt vs. counterion screening 

The parameter /o can be used to quantify the screening mechanism. Values close to 
unity (for a monovalent system) indicate a system where the amount of counterions in 
the system is dominant, while a value close to zero means that the salt ions dominate the 
screening mechanism. In order to focus on the screening mechanism (salt screening vs. 
counterion screening), one should keep k (or hiR) strictly constant, while varying only 
fo. Within our numerical approach, this is easily possible: The computer experiment 
consists of adding more and more salt (in terms of concentration), which enhances the 
screening, while at the same time the box size is increased, such that the counterions 
are more and more diluted and their contribution to the screening is reduced. This is 
done in such a way that the total amount of screening remains constant (see Eq. ISTll . 
Of course, the reduced charge and the reduced diffusion coefficients are kept constant 
as well. 

A common assumption is that only the screening length, but not the screening 
mechanism should contribute to the value of the electrophoretic mobility. Although 
Fig. [2] and previous studies [22] show that within the given accuracy the effect of fo on 
the mobility is at least weak, there is no fundamental reason why the mobility should 
be strictly independent of /q. 

In order to test this quantitatively, we have studied (i) a single colloidal sphere in 
a box, corresponding to a simple-cubic (sc) crystal, (ii) two spheres in the box, such 
that the resulting crystal is body-centered cubic (bcc), and (iii) four spheres arranged 
in such a way to construct a face-centered cubic (fee) crystal. The fixed parameters 
were R = 1, Z = 10, D = 1.5, d = 0.08, while /o was varied. Since we work at constant 
resolution, the resulting curves {fired vs. fo) in Fig. H] are smooth. Indeed, the plot 
nicely shows that the reduced mobility does depend on the screening mechanism in a 
non-trivial fashion, for all three types of lattice structures. However, the effect is only 
of the order of 5 to 10%, and hence was not observable previously. 

In order to assess the effect of the lattice resolution on this result, we analysed the 
sc case in some more detail by varying d as well. Due to computational limitations, 
this was however confined to a fairly narrow interval 0.06 < d < 0.09. Nevertheless, a 
reasonably reliable linear extrapolation to the continuum limit c? — )■ seems possible, see 
Fig. [5l where this is shown for the data point fo = 0.5. The result of this extrapolation, 
namely the reduced mobility as a function of fo in the limit c? — > 0, is presented in 
Fig. ini where cubic splines were used for interpolation. Even though the data may not 




Figure 5. Reduced electrophoretic mobility as a function of the sphere resolution d 
for a constant value of /o = 0.5. The red line shows a linear fit to the numerical data. 




Figure 7. Reduced electrophoretic mobility as a function of the diffusion constant. 
Both curves are generated for a system with R ~ 1, Z ~ 10, d = 0.07 and various 
values of /o as indicated by the legend. 
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Figure 8. Reduced mobility as a function of /o for one single sphere in a periodic box 
(sc lattice) and the same parameters as in Fig. |4] for different diffusion constants as 
indicated by the legend. This plot shows d extrapolated values. 



be fully reliable due to the smallness of the d interval, they nevertheless indicate fairly 
convincingly that the effect is more than just a discretisation artifact. 

Furthermore, it turns out that this non-trivial behaviour is even more interesting 
when studying the effect of the diffusion constants. The limiting behaviour for — )■ 
and for — )■ oo can be easily understood. For — )■ 0, only the convective part of the 
convection-diffusion equation remains: 



0. 



4) 



At least for a single sphere in infinite space it is easily shown that this enforces the 
trivial solution v^^^ = (and hence /i = 0): In a spherical coordinate system with 
origin at the center of the sphere, and polar angle relative to the driving field, it is clear 
that the radial component of v^^^ must vanish, due to the above condition, while the 
azimuthal component must be zero as well, due to symmetry. Therefore only the polar 
component remains. However, imposing the condition V ■ v^^'> = for that component 
yields a solution that is either singular (and hence forbidden) or zero. We hence find 
/i — for — > 0, and it is highly plausible that this holds in the general well. 
Conversely, the limit — t- oo means that the convective term can be ignored, and the 
problem becomes independent of Di. Hence the electrophoretic mobility saturates at 
some limiting value. 

These predictions are nicely confirmed in our calculations, see Fig. [71 The data 
shown there were calculated for R = 1, Z = 10 and d = 0.07, but for different screening 
mechanisms. Since the curves intersect, one sees that switching from salt-dominated 
to counterion-dominated screening may either enhance the mobility (this happens for 
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small diffusion coefficients) or reduce it (this is the behaviour at large D). In more 
detail, this behaviour is analysed in Fig. |8l 

4 ■4- Weakly charged colloids 

A very interesting phenomenon can be observed in the case of weakly charged colloids. 
Consider an uncharged spherical obstacle in a solution of negatively and positively 
charged ions. Applying a constant external electric field, electro-osmotic flow is 
generated by electric forces acting on the salt ions; positive charges move with the 
field direction, negative charged ions against it. Since the ions can not penetrate the 
solid sphere, the ion fluxes will be deflected by the particle. Thus, negative salt ions 
accumulate at one side of the sphere, while positive ions are depleted in that region. 
Since no electric forces act on the uncharged particle, the accumulation of positive ions 
at one side and negative ions at the opposite side must be symmetric. This accumulation 
effect leads to a polarisation of the system. Note that the induced dipole moment points 
in the "wrong" direction, i. e. anti-parallel to the driving field. This is interesting, 
because it is in contrast to the usual observation that for a charged sphere the induced 
dipole moment points in the "right" direction, i. e. parallel to the external field [21] . 

For an uncharged sphere in an infinite salt solution, the problem is dramatically 
simplified and amenable to an exact analytical solution; this was recently presented by 
Dhont and Kang [H]. In reduced units their result for the dipole moment is 

p = -27111^ E,,t . (85) 

Our numerical calculations nicely reproduce this prediction for a system with Z = 0, 
R = 1, E = 1. However, again we need to reach the limit of an infinite system, meaning 
$ — !■ 0, as well as the continuum limit c? — )■ 0. This double extrapolation is presented in 
Figs. |9] and [101 and our final numerical result is 

p~-6.23, (86) 

which deviates less than 1% from the expected value — 27r. 

In a second step, the reduced charge Z of the colloidal sphere was increased, while 
we kept R = 1, E = 1 and used D = 1.5. The results for the dipole moment are 
presented in Fig. [TTl and they will be discussed below. 

We would like to comment that we believe that at this point the advantages of 
our perturbative approach come to full effect. In a non-perturbative calculation, the 
induced dipole moment would be a very weak signal on top of the leading-order charge 
cloud, and this weak contribution may be fairly difficult to be disentangled from artifacts 
in the leading order (discretisation errors and roundoff errors which result in an artificial 
nonzero dipole moment). Conversely, in our calculation the leading-order and the first- 
order contribution are cleanly separated. We hence believe our method to be more 
accurate and stable than non-perturbative approaches. 

Returning to Fig. [HI we observe that the dipole moment increases with increasing 
charge. For a reduced charge of about Z = 2, depending on the volume fraction, the 




Figure 10. Dipole moment as a function of the resolution d, where the data points 
are the results of extrapolations to the L — > oo limit. 
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Figure 11. Dipole moment as a function of the reduced particle charge, for i? = 1, 
E — 1, and various volume fractions and resolutions as indicated. 




Figure 12. 2-dimensional cuts of the first-order concentration profile of the negative 
salt ions , for the parameters d = 0.08, = <i>~8.18-10~^ and an electric field 
of Ex ~ 1 acting in the x direction. The reduced charges of the sphere and the dipole 
moments of the system are (a) Z — 0.0, Px ~ —6.44, (b) Z = 1.8, px = 0.087 and (c) 
Z = 3.0, Px = 11.16. 
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Figure 14. 2-dimensional cuts of the first order charge density p^^\ The parameters 
are the same as in the previous images (see Fig. [T^ . 




Figure 15. Same data as Fig. [TU^b) but with rescaled colormap. 

sign of the dipole moment changes, i. e. the orientation of the charge cloud is reversed. 
In order to elucidate the phenomenon in some more detail, we present in Figs. [T21 [13] 
and [H] the first-order charge clouds as two-dimensional cuts in the xy plane (the field 
is oriented in x direction): Figure [T^ depicts the negative salt ions. Fig. [T3]the positive 
salt ions and Fig. [TUthe charge density. Most interesting are the figures in the case of a 
very small but already "normal" dipole moment: Here one sees that the orientation of 
the charge cloud near the colloid is still "anomalous" , but this is more than compensated 
from "normal" contributions further away. For the charge distribution, this cut is plotted 
again in Fig. [15] with a rescaled colormap for better visibility. All in all, this charge 
cloud reversal highlights that in electrophoresis both electrostatic and hydrodynamic 
effects are important, and that these may compete, resulting in a change of qualitative 
behaviour depending on conditions. 

The critical charge, i. e. the value of the colloid charge at which the dipole moment 
switches its sign, depends on the volume fraction. Keeping all other parameters fixed 
as before and the resolution at d = 0.07, we mapped out the functions p vs. Z and 
determined the critical value via spline interpolation. The dependence of this value on 
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Figure 16. Critical charge as function of the inverse box lengtli L ^ as obtained from 
the roots of spline functions. The line is a linear fit to the data points. 
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Figure 17. Dipole moment as function of the reduced charge for the same parameters 
as in Fig. [TT] Here, we choose d = 0.08 and $ ~ 8.18 • 10~^. The blue circles are the 
same data as before. For the red squares the convection velocity was set to ■w^^^ = 0. 
The dash-dotted lines are cubic splines. 



the linear system size L is shown in Fig. [161 extrapolation yields 

^ 0) = 1.11 ±0.02. (87) 

Finally, Fig. [17] demonstrates that the charge cloud reversal is fairly strongly 
affected by the value of the ion diffusion coefficient. To this end, we also studied the 
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case Di oo, which we reahzed computationally by just repeating our calculation 
with D = 1.5 but turning the convection term in the convection-diffusion equation 
off. Alternatively, one may therefore view this calculation as a study that elucidates 
the influence of convection. The result is clearly that convective transport helps in 
establishing the "normal" orientation. 

5. Concluding remarks 

In this paper we investigated a new numerical approach for the theoretical treatment of 
a charge-stabilised colloidal dispersion in an external electric field. The system, given 
by a solid charged sphere in an electrolyte solution, was treated on a Mean-Field level, 
resulting in a system of coupled nonlinear partial differential equations. Following the 
ideas of O'Brien and White [13], the nonlinearity is confined to the equilibrium Poisson- 
Boltzmann equation by application of a linearisation with respect to the external field. 
The iterative procedure in combination with the chosen specialised solvers turns out to 
be very efficient and only limited by the choice of the lattice spacing. While the demand 
on memory for the bulk methods increases linearly with the number of grid nodes, the 
surface integral solver for the Stokes equation requires a dense matrix connecting all 
surface nodes of the colloidal particles. Thus, the amount of memory needed for the 
storage of this matrix increases rapidly with the resolution of the sphere. Therefore this 
method is very efficient up to a certain value of the resolution; beyond that, alternative 
solutions must be developed. However, the iterative method has the advantage that it 
is designed as a modular solver and every module can be replaced via an alternative 
algorithm. One possible approach for going to higher resolutions is to replace the Stokes 
solver by a bulk method. For example, the time step of a lattice Boltzmann method could 
be adjusted such that it is identical with the time step of the convection-diffusion solver. 
Thus, the iterative method could be modified in a way that the Nernst-Planck and the 
Stokes equation are solved simultaneously. Furthermore, the lattice Boltzmann method 
would have the same degree of locality as the convection-diffusion equation solver, 
and hence parallelisation using a domain decomposition would be easily implemented. 
Nevertheless, the current single-processor implementation is quite efficient and reliable 
for a fairly satisfactory range of parameters. The numerical results agree reasonably 
well with various established results from the literature. 

Furthermore, all parameters in the method are controlled independently, which 
offered, e. g., the opportunity to study the dependence of the electrophoretic mobility 
on the diffusion coefficient of the surrounding ions. One of the most interesting results 
presented above is the conclusion that the screening mechanism has an effect on the 
electrophoretic mobility, i. e. the mobility varies by a few percent if the amount of salt 
is increased in the solution, while the screening length is kept constant. This shows 
clearly that the assumption that counterion-dominated systems may be mapped onto 
salt-dominated systems is only approximately true: If the accuracy goes beyond, say, 
5%, special care must be taken for the exact screening mechanism. Moreover, this 
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dependence is also affected qualitatively by the diffusion coefficient. 

Another very interesting application is the examination of weakly charged colloidal 
systems. If an electric field acts on an uncharged colloidal sphere in salt solution, ions 
move and are deflected at the surface of the solid particle, resulting in an "anomalous" 
dipole moment anti-parallel to the driving field. This "anomalous" dipole moment was 
recently addressed by Dhont and Kang [H] analytically. With our numerical method we 
were able to reproduce their result up to one percent difference. Increasing the colloid 
charge, we find a critical value at which the dipole moment changes its sign and the ion 
cloud reverses its orientation. 

All in all, the developed tool is computationally much cheaper than the raspberry 
MD/LB model [39], while having a somewhat broader range of applicability than the 
original work of O'Brien and White [13]. Clearly, it has limitations, as outlined in 
more detail in the Introduction. Although it will therefore not be able to study all 
electrokinetic phenomena in charge-stabilised colloidal dispersions — in particular, 
many-colloid systems where the particles continuously move with respect to one another, 
and a global rest frame does not exist, are out of reach for the present single-colloid 
version — we believe that it has already proven useful and is fairly likely to continue to 
do so. 
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